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ABSTRACT 



O ; 

<"*] Aims. The earliest phases of massive star formation in clusters are still poorly understood. Here, we test the hypothesis for high-mass star 
Oh formation proposed in our earlier paper (Peretto et al. 2006) stating that a massive, ultra-dense core may be currently forming at the center of 
^ , the collapsing NGC 2264-C protocluster via the gravitational coalescence of several intermediate-mass Class objects. 

Methods. In order to confirm the physical validity of this hypothesis, we carried out IRAM Plateau de Bure interferometer observations of 
NGC 2264-C and performed SPH numerical simulations of the collapse of a Jeans-unstable, prolate dense clump. A detailed comparison 
, between these hydrodynamic simulations and both our earlier IRAM 30m observations and the new interferometer observations is presented. 
^* ' Results. Our Plateau de Bure observations provide evidence for disk emission in three of the six Class 0-like objects identified earlier with the 
. ^ 30m in the NGC 2264-C clump. Furthermore, they reveal the presence of a new compact source (C-MM13) located only ~ 10000 AU away, 
i but separated by ~ 1.1 km.s 1 in (projected) velocity, from the most massive Class object (C-MM3) lying at the very center of NGC 2264-C. 
' Detailed comparison with our numerical SPH simulations supports the view that NGC 2264-C is an elongated cluster-forming clump in the 
process of collapsing and fragmenting along its long axis, leading to a strong dynamical interaction and possible protostar merger in the central 
region of the clump. The marked velocity difference observed between the two central objects C-MM3 and C-MM13, which can be reproduced 
in the simulations, is interpreted as an observational signature of this dynamical interaction. The present study also sets several quantitative 
constraints on the initial conditions of large-scale collapse in NGC 2264-C. Our hydrodynamic simulations indicate that the observed velocity 
pattern characterizes an early phase of protocluster collapse which survives for an only short period of time (i.e., < 1 x 10 5 yr). To provide 
a good match to the observations the simulations require an initial ratio of turbulent to gravitational energy of only ~ 5 %, which strongly 
suggests that the NGC 2264-C clump is structured primarily by gravity rather than turbulence. The required "cold" initial conditions may result 
from rapid compression by an external trigger. 

Conclusions. We speculate that NGC 2264-C is not an isolated case but may point to key features of the initial phases of high-mass star 
formation in protocluster. 

Key words, stars: formation - stars: circumstellar matter - stars: kinematics - ISM: clouds - ISM: kinematics and dynamics - Hydrodynamics 
- ISM: individual object: NGC 2264-C 



1. Introduction 



Most stars are believed to form in clusters (e.g. Lada & Lada 
2003) and high-mass stars may form exclusively in cluster- 
forming clouds. For a comprehensive understanding of clus- 
tered star formation, a good knowledge of the initial conditions 
and earliest phases of the process is crucial which, in practice, 
can only be inferred from detailed studies of deeply embedded 



protoclusters at (sub)millimeter wavelengths (e.g. Motte et al. 
1998, Andre 2002). 

Two main scenarios have been proposed to explain the 
formation of high-mass stars in clusters. In the first scenario, 
high-mass stars form essentially in the same way as low-mass 
stars, via an enhanced accretion-ejection phase. In the stan- 
dard model of low-mass star formation, the mass accretion rate 
is governed by the thermal sound speed and does not exceed 
~ 10~ 5 Mo.yr 1 (e.g. Shu 1977; Stahler et al. 2000) in cold 
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cores. To form high-mass stars by accretion, a significantly 
higher accretion rate ~ 1CT 3 M .yr~' is required to overcome 
the radiation pressure generated by the luminous central object 
(e.g. Wolfire & Cassinelli 1987). In order to solve this problem, 
McKee & Tan (2003) proposed a model in which high-mass 
star formation takes place in ultra-dense cores supported by tur- 
bulence within virialized cluster-forming clumps. This model 
produces high-mass accretion rates such as those required to 
form high-mass stars by accretion. 

In the second scenario, high-mass stars form by coales- 
cence of lower-mass stars in the dense inner core of a contract- 
ing protocluster (Bonnell et al. 1998). This scenario requires 
high stellar densities (i.e. ~ 10 8 stars.pc -3 ) in order to render 
the probability of stellar collisions high enough and allow stel- 
lar mergers to take place. It avoids the accretion problem of 
high-mass star formation by directly combining the masses of 
lower, intermediate-mass stars. However, no detailed model ex- 
ists yet to describe how this coalescence mechanism actually 
occurs. 

On the observational side, studying the earliest stages of 
high-mass star formation is particularly difficult due to the 
tight packing, deeply embedded nature, and relatively large 
distances of massive protoclusters. Based on IRAM 30m ob- 
servations of the massive cluster-forming clump NGC 2264- 
C (d~ 800 pc), Peretto, Andre, Belloche (2006 - hereafter 
PAB06) recently proposed a picture of high-mass star forma- 
tion combining features of the two above-mentioned scenar- 
ios. They showed that NGC 2264-C harbored a dozen Class 0- 
like objects (cf. Andre, Ward-Thompson, Barsony 2000) and 
was characterized by large-scale collapse motions (see also 
Williams & Garland 2002). They suggested that a massive, 
ultra-dense protostellar core was in the making in the cen- 
tral part of the NGC 2264-C clump as a result of the grav- 
itational merger of two or more lower-mass Class objects. 
The total mass inflow rate associated with the collapse of the 
clump toward the central protostellar core was estimated to be 
3 x 10~ 3 M .yr _I . PAB06 argued that the combination of large- 
scale collapse and protostellar mergers may be the key to pro- 
duce the conditions required for high-mass star formation in 
the center of NGC 2264-C. 

In this paper, we follow up on the detailed single-dish study 
of NGC 2264-C by PAB06 and present higher-resolution ob- 
servations of the same cluster-forming clump taken with the 
IRAM Plateau de Bure interferometer. We compare our obser- 
vations with a set of SPH hydrodynamic numerical simulations 
which attempt to specifically model NGC 2264-C. As the kine- 
matical and density patterns of NGC 2264-C appear to be rel- 
atively simple, comparison between millimeter observations of 
this region and numerical models offers a unique opportunity to 
make progress in our understanding of clustered star formation. 

Section 2 presents our PdBI observations. Section 3 de- 
scribes the dedicated hydrodynamic SPH simulations that we 
performed to model NGC 2264-C. We compare the observa- 
tions with the numerical simulations in Sect. 4 and draw several 
conclusions in Sect. 5. 



2. Interferometer observations of NGC 2264-C 

2.1. Observations 

We performed 3.2 mm and 1.2 mm observations of the cen- 
tral part of NGC 2264-C with the IRAM Plateau de Bure in- 
terferometer (PdBI) in December 2003 and April 2004. We 
used the C and D configurations with 6 antennas. We used 
both 1 mm and 3 mm receivers with 244.935620 GHz (A = 
1.2 mm) and 93.176258 GHz (A = 3.2 mm) as central rest 
frequencies. We observed at four positions which were cho- 
sen so as to obtain a fully sampled mosaic at 3.2 mm (pri- 
mary beam FWHM ~ 54") and to encompass the millimeter 
sources C-MM1, C-MM2, C-MM3, C-MM4, CMM5, and C- 
MM9 (see Fig. 2b of PAB06) identified by PAB06 with the 
IRAM 30m telescope. Because the corresponding 1.2 mm mo- 
saic is undersampled (primary beam FWHM ~ 20"), only two 
of these sources (C-MM3 and C-MM4) were effectively im- 
aged at 1.2mm. We obtained a 3.2 mm dust continuum mo- 
saic and two separate 1.2 mm continuum maps, as well as a 
N 2 H + (l-0) mosaic. The spectral resolution for the N 2 H + (l-0) 
data was 20 kHz, which corresponds to a velocity resolution of 
0.06 km.s 1 at 93.2 GHz. The sources used for the bandpass, 
amplitude and phase calibrations were 0420-014, 0528+134, 
0736+017, CRL618, 0923+392, and 3C273 (only the first three 
were used for the second run in April). We calibrated the 
data and produced images using the CLIC and MAPPING 
softwares (Lucas 1999, Guilloteau et al. 2002), part of the 
GILDAS packag43 (Pety 2005). The deconvolution was per- 
formed using the natural weighting option of the Clark (1980) 
CLEAN algorithm (Guilloteau 2001). The final synthesized 
beamwas5.3"x3.8" (HPBW) with PA.=+63° at 3.2mm, and 
3.1"xl.5" (HPBW) withP.A.=+74° at 1.2 mm. We also com- 
bined our PdBI N2H + (l-0) observations with the single-dish 
N 2 H + data cube of PAB06 in order to recover short-spacing 
information. The resulting synthesized beam of the combined 
N 2 H + (l-0) mosaic is 6.1"x4.0" (HPBW) with PA.=+65°. 

2.2. Dust continuum results: Evidence for disk 
emission 

Our PdBI 3.2 mm dust continuum mosaic is shown in Fig. Q] 
It reveals only pointlike sources, since most of the extended 
emission seen in the single-dish 1 .2 mm dust continuum map 
(cf Fig. 2 of PAB06) was filtered out by the interferometer. The 
final rms noise level was <x ~ 0.8 mJy/beam. We extracted mil- 
limeter sources from this map using the Gaussclump algorithm 
(Stutzki & Giisten 1990). We detected seven peaks lying above 
5<t. All of these peaks were previously detected with the 30m 
telescope, except one object, here called C-MM13 (cf. Fig. 1), 
which is a new detection. Several other peaks lie between 3<x 
and 5cr, but, by lack of confidence in these marginal detections, 
we did not consider them. The present 3.2 mm continuum map 
confirms and improves the positions of the compact millimeter 
continuum sources detected by PAB06. Interestingly our PdBI 
1 .2 mm continuum observations of C-MM3 and C-MM4 (not 

1 See http : //www . iram . fr/IRAMFR/GILDAS for more informa- 
tion about the gildas softwares. 
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Table 1. Measured source properties 



Source 


Coordinates 


Undec.FWHM 


PA. 


01.2 

J peak 


c 30/n 
peak 


peak 


03-2 
peak 


c3.2 




(«2000 £2000) 


(arcsecond) 


(deg) 


(mjy/beam) 


(mJy/beam) 


(mJy/beam) 


(mJy/beam) 


(mJy) 


[1] 


[2] 


[3] 


[4] 


[5] 


[6] 


[V] 


[8] 


[9] 


C-MM1 


06:41:17.95 +09:29:03 


5.8x4.1 


63 




255 


6 


18 


21 


C-MM2 


06:41:15.15 +09:29:10 


5.6x4.3 


50 




183 


4 


7 


8 


C-MM3 


06:41:12.30 +09:29:12 


5.7x4.4 


58 


224 


573 


14 


37 


45 


C-MM4 


06:41:09.95 +09:29:22 


7.0x5.0 


87 


77 


426 


10 


14 


24 


C-MM5 


06:41:10.15 +09:29:36 


6.2x4.5 


83 




261 


6 


5 


7 


C-MM9 


06:41:15.30 +09:29:33 


7.3x4.3 


51 




94 


2 


6 


9 


C-MM13 


06:41:11.45 +09:29:17 


7.5x6.1 


87 








5 


11 



[1]: The C-MM numbers are the same as in PAB06. The new source is labelled C-MM13. 

[2]: J2000 source coordinates, accurate to better than 1", derived from a Gaussian fit to the PdBI 3.2mm dust continuum map. 

[3]: Undeconvolved FWHM sizes derived from fitting an elliptical Gaussian to the PdBI 3.2mm dust continuum map. 

[4]: Position angle (from North to East) of the major axis of the fitted Gaussian ellipse from the 3.2 mm dust continuum map. 

[5]: PdBI 1.2mm peak flux density at the best-fit source position (HPBW=3.1"xl.5") 

[6]: 30m 1.2mm peak flux density at the source position (HPBW=11"; from PAB06) 

[7]: 3.2mm peak flux density expected at PdBI angular resolution estimated from col. [6] (HPBW=4.5") 

[8]: PdBI 3.2mm peak flux density at the best-fit source position (HPBW=5.3"x3.8") 

[9]: PdBI 3.2mm integrated flux density inside the fitted Gaussian ellipse 



shown here) do not reveal any further sub-fragmentation. The 
source properties as derived from our PdBI and 30m observa- 
tions are summarized in Table Q] 

The millimeter dust continuum emission from a (Class C 
or Class I) protostar a priori originates from two components: 
an extended envelope (a few thousand AU in size) and a disk 
(up to a few hundred AU). At a distance of 800 pc, the lin- 
ear resolution of the 30m telescope is ~ 9000 AU (HPBW) at 
A = 1.2 mm. On this spatial scale, the 1.2mm continuum emis- 
sion observed toward a young protostar (TableQ]col.[6]) is ex- 
pected to be dominated by the envelope rather than by the disk 
(e.g. Andre, Ward-Thompson, Barsony 2000; Looney, Mundy. 
Welch 2000). Conversely, we expect the disk component tc 
dominate on compact, interferometric scales (e.g. Terebey et 
al. 1993). Assuming an isothermal, centrally-condensed enve- 
lope, i.e. with a density p oc r -2 , the flux density is expected 
to scale linearly with beam size: S v oc 6. We thus expect the 
envelope contribution to the PdBI flux density at 1 .2 mm to be 
given by: 

„exp:\.2 _ „ 3 ()m v ( HPBW Bure \ m 

If we assume the Rayleigh-Jeans regime and adopt a dust 
opacity scaling as k v oc V s (e.g. Hildebrand 1983), then we can 
also estimate the expected contribution of the envelope to the 
PdBI flux density at 3.2 mm: 

/ 1 7\^+2 

„exp:3.2 _ „exp:1.2 I i n ^ 
peak ~ peak I 32 I 

In order to estimate a lower limit to the disk component, 
we choose /? = 1 which maximizes the contribution of the en- 
velope. A value of /? = 1 .5 is likely to be more representative of 
protostellar cores/envelopes (e.g. Ossenkopf & Henning 1994) 
and would yield a lower estimate for the expected envelope 
contribution. The expected envelope flux densities are listed in 



3.2mm dust continuum (PdBI) 




50 o -50 

Aa (arcsec) 



Fig. 1. 3.2 mm dust continuum mosaic of the central part of 
NGC 2264-C obtained with the PdBI. The (0", 0") position 
corresponds to (a = 06*13 m 00*; 6 = 09°29'10") in J2000 co- 
ordinates. The pre/protostellar sources detected in this map are 
marked with white crosses and labelled C-MM, keeping the 
same numbering scheme as in PAB06. The two black crosses 
mark the positions of C-MM 10 and C-MM 12 (cf. PAB06), 
which are located at the very edge of our mosaic. The black 
open star shows the position of the luminous IRAS source 
IRS1. The rms of the mosaic is <x =s 0.8 mJy/beam. The 
contours are 2.4 mJy/beam (i.e. 3(f), 4 mJy/beam (i.e. 5<t), 
5, 6, 8 mJy/beam and then go from 10 to 35 mJy/beam by 
5 mJy/beam. 



Table[T]col.[7] for each detected source. It can be seen that they 
are a factor of ~ 2 - 3 lower than the observed flux densities 
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(TableQ]col.[8]) for C-MM1, C-MM3 and C-MM9. The excess 
flux density observed on small spatial scales can be attributed 
to unresolved disk emission (e.g. Terebey et a. 1993). Our re- 
sults thus suggest the presence of a disk in C-MM1, C-MM3, 
C-MM9 and confirm the protostellar nature of these candidate 
Class sources. Given the uncertainties on the dust emissivity 
index (3 (e.g. Dent et al. 1998), we cannot conlude on the pres- 
ence or absence of a disk in C-MM2 and C-MM4. Finally, it is 
very likely that C-MM5 does not have a disk since its observed 
3.2 mm flux density is consistent with pure envelope emission. 

For the three sources showing evidence of disk emission, 
i.e. C-MM1, C-MM3 and C-MM9, we can estimate both the 
disk and envelope masses as follows. First, we estimate the 
flux arising from the disk by subtracting the expected envelope 
peak flux density given in Table Q] col. [7] from the observed 
peak flux density of col. [8] (here, we assume that our 3.2 mm 
PdBI observations do not spatially resolve the disk given the 
distance of NGC 2264). Then, the flux arising from the enve- 
lope is considered to be given by the integrated flux of col. [9] 
minus the disk contribution. These flux estimates are listed in 
Table [2] For the envelope mass estimates, we follow PAB06 
and assume a dust temperature T,/ = 15 K, J3 = 1.5, as well 
as a dust opacity K\.2mm = 0.005 cm 2 .g _1 corresponding to 
K3.2mm = 1.3 x 10~ 3 cm 2 .g _1 . Concerning the disk mass esti- 
mates, we assume a dust temperature of T^ = 50 K and a dust 
opacity K\2mm - 0.02 cm 2 .g _1 (Beckwith et al. 1990) corre- 
sponding to Kj,.2mm = 5.2 x 10~ 3 cm 2 .g _I (see Table |2). The 
dust temperature is higher for the disk because it is supposed to 
be warmer, closer to the star, while the dust opacity is slightly 
different from the one adopted in PAB06 because of the en- 
hanced dust emissivity expected in the dense central parts of 
protostellar disks (e.g. Ossenkopf & Henning 1994). We cau- 
tion that the disk masses calculated in this way are only rough 
estimates. A more proper analysis of the density structure of 
these sources, especially C-MM3, through sub-arcsecond mil- 
limeter observations would be of great interest. By deconvolv- 
ing the FWHM sizes of Table Q] from the synthesized beam, 
we derive the geometrical mean diameter of each source. For 
each object, we can then estimate the mean column density and 
mean volume density of the envelope component. All of these 
derived source parameters are listed in Table [3] It can be seen 
that C-MM3 and C-MM4 have the most massive envelopes 
by far, with M em > 15 M Q in both cases. The densest enve- 
lope/core is associated with the central source C-MM3, which 
reaches a mean volume density > 1x10 s cm -3 on a ~ 3200 AU 
(2XFWHM) scale. 

The new source identified with the interferometer, C- 
MM13, could not be separated from C-MM3 at the angular 
resolution of the IRAM 30m telescope. The projected dis- 
tance between C-MM3 and C-MM13 is only ~ 10000 AU. 
Surprisingly, C-MM13 is the only new compact millimeter 
continuum source detected with the PdBI above 5<x. The fact 
that there is almost no sub-fragmentation despite a factor 2- 
4 improvement in angular resolution between the partially re- 
solved 30m sources and their PdBI counterparts suggests that 
most of the compact dust continuum sources detected at the 
30m represent individual Class objects rather than small 
groups of protostellar cores. However, we stress that our PdBI 



N 2 H + (1-0) Integrated intensity map (PdBI) 




beam 



50 -50 

(arcsec) 

Fig. 2. N2H + (l-0) integrated intensity map of the central part 
of NGC 2264-C obtained with the PdBI. The (0", 0") posi- 
tion corresponds to (a = 06*13 m 00 s ; 6 = 09°29'10") in J2000 
coordinates. The crosses with labels show the positions of the 
pre/protostellar cores detected in our PdBI dust continuum mo- 
saic (Fig. [TJ, while the black cross without label mark the po- 
sition of C-MM10 (cf. PAB06). The white open star shows the 
position of the IRAS source IRS 1 . The contours go from 1 to 5 
Jy/beam.km/s by 1 Jy/beam.km/s. 

Table 2. Estimated PdBI 3.2 mm flux densities of the disk and 
envelope components for the three objects showing evidence of 
disk emission 



Source 


c a 
(mJy) 


peak — 

(mJy/beam) 


Cenvc 

(mJy) 


(M Q ) 


C-MM1 


12 


6 


9 


0.6 


C-MM3 


23 


14 


22 


1.1 


C-MM9 


4 


2 


5 


0.2 



" 3.2 mm flux density of the disk component estimated by subtract- 
ing the peak flux density expected for the envelope (col. [7] of TableQ} 
from the observed peak flux density (col. [8] of TableQ} 

b Estimated 3.2 mm peak flux density of the envelope component at 
the PdBI resolution (cf. col.[7] of TableQ} 

c 3.2 mm integrated flux density of the envelope component esti- 
mated by subtracting the disk contribution given in col. [2] from the to- 
tal integrated flux density measured with PdBI (cf. col. [9] of TableQ} 

d Disk mass seen with the PdBI and estimated from SdM (with Tj = 
50 K and k = 5.2 x 10~ 3 cm 2 .g~'). Typical uncertainty is a factor > 2 
(on either side) due to uncertain dust opacity and dust temperature. 



observations could only detect dust continuum sources with 
3.2 mm peak flux densities larger than 4 mJy/beam (i.e. 5cr). 
Assuming the same temperature and dust properties as above 
for the envelope mass estimates, this corresponds to a mass 
detection threshold ~ 3 M . Therefore, we could not have 
detected low-mass pre/protostellar cores possibly lying in the 
vicinity of the main sources listed in TableQ] This point is fur- 
ther discussed in the next section. 

2.3. N 2 H + (1-0) results 

2.3.1 . PdBI observations only 

Our PdBI N2H + (l-0) integrated intensity map is shown in 
Fig. [2] As for the dust continuum mosaic, most of the extended 
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Table 3. Derived source parameters 



Source 


FWH1VF 




M em £ 










(AU)~ 


(10 23 cm" 2 ) 


(Mo) 


(Mo) 


(cm" 3 ) 


(km.s- 1 ) 


C-MM1 


1400 


7 


6.2 


13.1 


8.3xl0 7 




C-MM2 


1400 


9 


5.5 


16.0 


7.3xl0 7 


6.2 


C-MM3 


1600 


17 


15.2 


40.9 


1.4xl0 8 


7.1 


C-MM4 


3000 


17 


16.5 


35.1 


2.2xl0 7 


8.9 


C-MM5 


2200 


6 


4.8 


18.4 


1.6xl0 7 




C-MM9 


2400 


2 


3.5 


6.6 


9.3xl0 6 




C-MM13 


4000 


6 


7.6 




4.3xl0 6 


8.2 



" Geometrical mean of the deconvolved, major and minor FWHM diameters measured on the PdBI 3.2mm continuum mosaic 

* Column density of the envelope seen with the PdBI and estimated from S 3 2 n;t for the sources without a disk (col. [8] of TablefTJ and from 

Speak f° r tne sources with a disk (col. [3] of Tabled. The used dust properties are = 15 K and k = 1.3 x 10~ 3 cm .g . Typical uncertainty is 

a factor > 2 (on either side) due to uncertain dust opacity and dust temperature. 
Envelope mass seen with the PdBI and estimated from S™ 1 ' (Tabled with T rf = 15 K and k = 1.3 x 10 3 cm 2 .g-'). The uncertainty is the 

same as for the column density. 

d Core mass seen with the 30m telescope (from PAB06). This mass is larger than M em , because, contrary to the PdBI, the 30m does not filter 
out most of the large scale emission 

e Envelope volume density estimated in a radius equal to the FWHM. Estimated in a radius twice smaller, the density increases by a factor 
of 4. The uncertainty is the same as for the column density 

' LSR velocity estimated from the combined N2H + (l-0) spectra 



emission is filtered out. On this map, it can be seen that C- 
MM5 and possibly C-MM12 are not closely associated with 
a N 2 H + (l-0) peak. We also note that the strongest N 2 H + (l-0) 
peak is associated with the new C-MM13 object which is one of 
the weakest dust continuum sources. The morphological differ- 
ences between the dust continuum sources and their N2H + (l-0) 
counterparts may possibly reflect differences in chemical evo- 
lutionary stage (cf. Aikawa et al. 2005). It is also noteworthy 
that several N2H + (l-0) peaks are not associated with any of the 
bona fide dust continuum sources listed in Table[3] On the other 
hand, most of these N2H + (l-0) peaks have faint dust continuum 
counterparts with flux levels between 3<x and 5<x in the 3.2 mm 
continuum map (compare Fig.[TJand Fig. 

The N2H + (l-0) spectra observed at PdBI (prior to combi- 
nation with 30m data) provide an estimate of the velocity dis- 
persion within the sources on a ~ 3500 AU (FWHM) spatial 
scale. The mean line-of-sight velocity dispersion is found to 
be < (Ti os >- 0.34 km.s-'. Assuming a kinetic temperature 
T* = 15 K, the non-thermal contribution to this velocity dis- 
persion is < >^ 0.33 km.s-'. Comparing this value with 
the isothermal sound speed, c s 0.23 km.s" 1 , we conclude that 
the sources are still marginally dominated by non-thermal mo- 
tions (due to, e.g, turbulence, collapse or, outflow) on scales of 
a few thousand AUs. 

2.3.2. Combined PdBI and 30m observations 

As mentioned in § 2.1, we added short-spacing 30m informa- 
tion to our interferometric data in order to obtain a N2H + map 
sensitive to a wide range of angular scales from ~ 4" up to ~ 4'. 
This combination was performed using the MAPPING package 
developed by IRAM (Guilloteau et al. 2002). The combined 
PdBI/30m map of N2H + (l-0) integrated intensity is shown in 
FigEl As expected, more extended emission is present in the 



N 2 H + (1-0) integrated intensity map (30m + PdBl) 




Aa (arcsec) 

Fig. 3. N2H + (l-0) integrated intensity map of the central part 
of NGC 2264-C resulting from the combination of our 30m 
and PdBI data. The (0", 0") position corresponds to {a = 
06 A 13"W; 5 = 09°29'10") in J2000 coordinates. The crosses 
mark the positions of the pre/protostellar 3.2 mm dust contin- 
uum sources detected in Fig.Q] For the sake of readability, only 
a few sources are labeled. The white open star symbol shows 
the position of the IRAS source IRS 1 . The contours go from 3 
to 13 Jy/beam.km/s by 2 Jy/beam.km/s. 



combined mosaic, but the compact sources detected in the PdBI 
map are still clearly visible. 

In order to constrain the kinematical pattern of these 
sources within the NGC 2264-C clump, we constructed a 
position-velocity (PV) diagram along an East- West axis go- 
ing through the four central sources, C-MM2, C-MM3, CMM4, 
and C-MM13 (see Fig. 3J. This PV diagram shows an overall 
velocity gradient of 8.4 km.s" 1 .pc" 1 from East to West between 
C-MM2 and C-MM4. The LSR velocities of each of the four 
sources C-MM2, C-MM3, C-MM4, and C-MM13 are listed 
in Table [3] Figure |4] helps to clarify the origin of the veloc- 
ity discontinuity identified by PAB06 with the 30m telescope 
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Fig. 4. Position-velocity diagram derived from the combined 
(PdBI + 30m) N 2 H + (101-012) data cube by taking a cut along 
the East- West axis going through C-MM2, C-MM3, C-MM4 
and C-MM13. The positions of these sources are plotted on the 
position axis. The white contours go from 0.5 to 1.1 Jy/beam 
by 0.2 Jy/beam. 



in the center of NGC 2264-C (see Fig. 6 of PAB06). At the 
30m resolution, the N2H + (101-012) spectrum observed toward 
the central source C-MM3 was double-peaked. The higher res- 
olution of the PdBI interferometer now allows us to identify 
a distinct component, C-MM13, separated by 13" in position 
and ~ 1.1 km.s 1 in velocity from C-MM3. With the 30m tele- 
scope, C-MM13 could not be separated from C-MM3. We also 
observe in the western part of the PV diagram (at an offset of 
-70" and velocity of 7.4 km.s~') a velocity feature which is 
associated with the strong N2H + (l-0) peak lying in this part 
of the clump (cf Fig. [3). This velocity feature clearly departs 
from the rest of the diagram. We attribute it to a peculiar veloc- 
ity field around the luminous young star IRS 1 , whose wind has 
likely perturbed the ambient velocity field and triggered star 
formation in the immediate vicinity (~ 10" in radius around 
IRS 1), as suggested by the observations of Nakano et al. (2003) 
and Schreyer et al. (2003). 

To summarize, our interferometric observations confirm the 
(Class 0) protostellar nature of C-MM1, C-MM3, C-MM9, and 
set new constraints on the kinematics of NGC 2264-C. The 
PdBI observations can help us to confirm or refute the scenario 
proposed by PAB06 of an axial collapse of NGC 2264-C along 
its long axis leading to the merging of dense cores in the center. 
In the context of this scenario, the protostellar nature of the mil- 
limeter continuum cores sets strong timescale constraints: the 
individual collapse of the cores must occur on a significantly 
shorter timescale than the larger-scale collapse of the clump as 
a whole. The presence of two central sources, C-MM3 and C- 
MM13, adjacent to one another (i.e. 10000 AU) and with a ve- 
locity difference of ~ 1.1 km. s _1 , must also be accounted for. In 
the next two sections, we attempt to match these observational 
constraints with hydrodynamic simulations. 



3. SPH numerical simulations 

3.1. Numerical method and initial conditions 

In order to test the physical validity of the scenario pro- 
posed by PAB06 (see also § 1), we performed Smooth Particle 
Hydrodynamics (SPH) simulations (Monaghan 1992, Bate et 
al. 2003) using the DRAGON SPH code from the Cardiff insti- 
tute (see Goodwin et al. 2004). We simulated the time evolution 
of an isothermal (T* = 20 K), Jeans-unstable elongated clump 
of mass 1000 M , comparable to the estimated total mass of 
NGC 2264-C (~ 1600 M Q - PAB06). The model clump was 
initially ellipsoidal (finite boundary conditions) with an aspect 
ratio of 2. The initial density profile was such as: 

" H2 (1 + (r/r ) 2 + (z/2r ) 2 ) W 

corresponding to a flat inner (r < ro) region and a oc r~ 2 
outer (r > ro) region. The total mass of the flat inner region was 
~ 200 M . Highly concentrated clouds (with, e.g., n# 2 oc r~ 2 ) 
are known to hardly fragment during collapse (cf. Myhill & 
Kaula 1992; Whitworth et al. 1996), while uniform or mod- 
erately concentrated clouds (with n#, oc r _1 or flatter) typi- 
cally fragment into as many Jeans-mass fragments as they con- 
tain initially (cf. Burkert et al. 1997). The expected number of 
fragments produced by the collapse of our model clump thus 
corresponds to the number of Jeans masses estimated in the 
flat inner core. At the mean density calculated within ro (i.e. 
n//, ~ 1 x 10 3 cm 3 ) and for Tj = 20 K, the Jeans mass is 
My ^ 20 M (see Bonnell et al. 1996 for a precise defini- 
tion of M/), which yields a Jeans mass number Nj ~ 10 for 
the flat inner core. In these simulations we also included tur- 
bulent fluctuations. Since the exact nature and properties of 
interstellar turbulence are not fully understood yet, we con- 
sidered two types of energy spectrum: 1) a spectrum scaling 
as Kolmogorov turbulence, i.e., E(k) oc k~ 5 ^, and 2) a white 
spectrum, i.e, E(k) oc k°. The phases of these turbulence fluc- 
tuations were chosen randomly. Initially, three energy compo- 
nents controlled the evolution of the model filament: the ther- 
mal energy, 7~,/,, which remained constant in time throughout 
the simulations (i.e., isothermal assumption); the gravitational 
energy, <W, whose initial value depended on the clump den- 
sity profile and thus on n c and ro (larger values of n c and/or ro 
correspond to lower gravitational energy; cf. Eq.fO); and the 
turbulent energy, T tur b- In all the simulations shown in this arti- 
cle, n c and ro have the same value, namely ndHi) = 2000 cm -3 
and ro = 0.7 pc, corresponding to an initial thermal to gravita- 
tional energy ratio op, ~ 8%. The initial value of the ratio of 
turbulent to gravitational energy, oP , , was varied from 0% to 
50%. We have also explored other initial conditions. When n c 
is too large (n c > 5000 cm -3 ), we find that too many fragments 
form, when n c is too small (n c < 500 cm 3 ), we find that only 
one central fragment will generally form. We have also varied 
the initial aspect ratio and conclude that if it is too close to one, 
then the cloud is not sufficiently filamentary whereas it is too 
much filamentary if it is initially too large. Finally, we have also 
explored the possibility that NGC 2264-C could be the result 
of a collision between two preexisting clouds rather than a col- 
lapsing elongated clump. However, it was not possible to repro- 
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Fig. 5. Time evolution of a simulation with an initial turbulent to gravitational energy ratio a" turh = 5%. The first column displays 
the velocity field of the particles taken along the long axis (z-axis) of the model filament, as traced by individual SPH particles. 
The second column similarly displays the evolution of the density cut along the z-axis of the model filament. The third column 
shows synthetic column density maps in the (z,x) plane. The reference time, labelled t//, is taken to be one global free-fall time 
of the initial clump after the start of the simulation, i.e., t// = 9.5 x 10 5 yr. The first row is taken at a time step tyy - 4 x 10 s yr, 
while the second, third, and fourth rows are for time steps t// - 2 X 10 s yr, t// - 1.0 x 10 5 yr, and t// - 0.4 x 10 s yr, respectively. 



duce the various features of this cloud within the scope of this 
scenario. The most important disagreements are i) a collision 
quickly tends to create a sheet rather than a filamentary object; 
ii) we find it very difficult to produce a well defined third object 
like C-MM3 by interaction of two colliding clouds; iii) a sim- 
ple collision is unable to create a series of young condensations 



spread over the long axis of the clump at, e.g., positions com- 
parable to C-MM1 and C-MM5. All our simulations were per- 
formed with a total of 5 million SPH particles. When the local 
density exceeded n# 2 = 1.3 x 10 8 cm -3 , standard SPH particles 
were replaced by sink particles. The radius, r sir & = 500 AU, 
of the sink particles defines the highest resolution reached by 
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our simulation. All particles falling within r „„^. of a sink parti- 
cle and being bound to it were removed from the simulations, 
and their mass, linear and angular momentum were added to 
the corresponding sink particle values. Using sink particles al- 
lowed us to avoid artificial fragmentation (Truelove et al. 1997, 
Bate & Burkert 1997). The relatively low density threshold at 
which sink particles were introduced implies that we could not 
model advanced phases of star/cluster formation but only the 
first stages of clump fragmentation. Indeed, the limited numer- 
ical resolution of our simulations prevented us from describing 
small spatial scale processes such as disk formation. 

3.2. General Pattern 

Figure [5] displays the density and velocity fields along the z- 
axis (i.e., long axis), as well as the column density maps in 
the (z,x) plane, at four time steps for a model filament with an 
initial level of turbulence, or ' . = 5%. The reference time was 
chosen to be at one global free-fall time, tff, after the start of 
the simulation. Given the initial central density of the model 
clump, this corresponds to t// = 9.5 x 10 5 yr. The four time 
steps shown in Fig. [5] were taken at tyy - 4 x 10 s yr, t//-2x 
10 5 yr, tff - 1.0 x 10 5 yr, and t ff - 0.4 x 10 5 yr, respectively. 

We can describe the evolution of the model clump as fol- 
lows (see also Bonnell et al. 1996; Inutsuka & Miyama 1997). 
Since the ellipsoidal clump initially contains several thermal 
Jeans masses (i.e. Ny ~ 10, cf. Sect. 3.1) and has a shorter 
dynamical timescale perpendicular to its major axis, it first col- 
lapses along its minor axis when seen projected onto the plane 
of the sky (see first and second panels of Fig. [5j, amplifying 
the initial anisotropy and leading to the formation of a very 
elongated and filamentary structure (cf. Lin et al. 1965). This 
fast contraction proceeds until thermal and turbulent pressure 
gradients can stop the collapse along the minor axis, ensuring 
an approximate hydrostatic equilibrium in the (xy) plane (cf 
Bonnell et al. 1996). 

Since the dynamical timescale is longer along the major 
axis, the clump keeps collapsing along its long axis, i.e. the z- 
direction, after a transverse equilibrium has been established. 
The velocity field is initially nearly homologous (i.e., V z oc z- 
see first and second panels of Fig. |5]) but becomes more and 
more complex as the filament fragments into several cores, 
each of them collapsing individually as can be seen on the third 
and fourth panels of Fig. [5] The individual collapse of the cores 
leads to the formation of local protostellar accretion shocks and 
associated protostars. Moreover, due to the global collapse of 
the model clump toward its center, we can also see the forma- 
tion of a central shock which separates the eastern (z > 0) and 
western (z < 0) sides of the clump. Altogether this dynamical 
evolution produces a complex density and velocity pattern. 

When oft. b = 50%, the velocity field (not displayed here 
for conciseness) is much less organized. More shocks develop 
which lead to enhanced clump fragmentation through the pro- 
cess now widely referred to as "turbulent fragmentation" in the 
literature (e.g. Padoan & Nordlund 2002, Klessen et al. 2005). 
The number of shocks is larger and the model clump becomes 
more substructured as the initial level of turbulence increases 



(see also Jappsen & Klessen 2004). At the other extreme, the 
case with no initial turbulence at all, eft, = 0%, does not 
yield any fragmentation due to the lack of initial fragmentation 
seeds, which is not realistic. 

Therefore, based on the results of our simulations, several 
general conclusions can be drawn. In particular, the initial level 
of turbulence appears to play a key role for the global aspect 
of the model clump. The higher the turbulence, the more dis- 
persed and less filamentary the model clump is. For low levels 
of initial turbulence, i.e., low values of a? ,, the structure and 
kinematics of the clump are dominated by gravity, while for 
high levels of turbulence the clump is primarily structured by 
turbulence. 

4. Detailed comparison between the observations 
and the SPH simulations 

We performed a wide set of SPH simulations with different ini- 
tial parameters, e.g., different values of the initial level of tur- 
bulence, cft rb , and of the thermal to gravitational energy ratio, 
oft h . We did not find it necessary to use different initial turbu- 
lent velocity fields since the time evolution of the model clump 
depends only weakly on this. When calculating synthetic ob- 
servations, we varied the inclination angle of the long axis of 
the model clump with respect to the line of sight. An inclina- 
tion angle of 45 degrees was adopted to produce the synthetic 
images and diagrams shown in Fig. [6] and Fig. [7] For each set 
of simulations, we seeked a particular time step at which the 
synthetic data, when convolved to the resolution of the obser- 
vations, best matched the existing 30m and PdBI constraints. In 
the next subsection, we present our "best-fit" simulations and 
discuss the consequences of changing the best-fit parameters. 

4.1. Overall morphology: A fragmented filament 

The first important feature which must be reproduced by the 
simulations is the elongated shape of NGC 2264-C and the 
presence of several protostellar sources, lining up along the 
long (East- West) axis of the clump. 

Figure |6]compares the observed column density maps (first 
row) with synthetic maps obtained from simulations with four 
different initial levels of turbulence, i.e., eft . — 1% (second 
row), eft , =5% (third row), eft , = 20% (fourth row), and 
eft rb = 50% (fifth row), all assumed to be "observed" with a 
viewing angle of 45 degrees. The first column of Fig. [6] dis- 
plays the synthetic data convolved to the 30m resolution while 
the second column displays the data convolved with the PdB 
interferometer beam. When convolving the simulated data to 
the PdBI resolution, we included the effect of interferometric 
filtering so as to allow more direct comparison with the ob- 
servations. For this purpose, we used the UVJVIODEL task of 
the GILDAS package. This task generated a set of visibilities 
in the UV plane by calculating the values of the Fourier trans- 
form of the simulated input image at the observed UV base- 
lines. The rest of the points in the UV plane was set to zero. 
This method had the consequence of filtering out all extended 
structures present in the numerical simulations. 
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The four simulations shown in Fig. [6] are compared when 
the synthetic column density maps convolved to the 30m res- 
olution match the observed map best, except for the case with 
a< turb ~ -^% (cf. fifth row) which does not exhibit a filamen- 
tary shape at any time step in the simulation. The correspond- 
ing time steps all lie in the range between t// - 5 x 10 5 yr 
and tff - 1.5 x 10 5 yr. The total mass accreted onto sink par- 
ticles at these time steps ranges from 0% to 0.2% of the initial 
clump mass. As already mentioned, we are thus looking at the 
very first stages of the formation of a protocluster, i.e., when 
the first pre-/proto-stellar cores with typical mean volume den- 
sities ~ 10 5 cirT 3 appear. It can also be seen in Fig|6]that the 
filamentary, elongated appearance of the NGC 2264-C clump 
cannot be reproduced when the initial level of turbulence is too 
high in the model. Based on this argument, we conclude that 
or?, has a maximum value of 20%, although the 5% model al- 
ready provides a better match to the observations than the 20% 
model. Figure[6]also shows that the large-scale morphology of 
the clump observed at the resolution of the 30m telescope pro- 
vides the strongest discriminator between different values of 
a 

Our "best-fit" simulation is shown in the third row of Fig. [6] 
It corresponds to a flat energy spectrum, E(k) cc k°, and an 
initial value of aft , =5%, which is much lower compared 
to other numerical SPH studies of cloud fragmentation (e.g. 
Bate et al. 2003). Note that the case of Kolmogorov-like turbu- 
lence leads to results which are broadly similar to that shown 
in Fig. [6] except that the shape of the filament is more irregular 
and too distorted to match the observations well. We therefore 
restrict our attention to the E(k) oc k° case. Although this en- 
ergy spectrum differs from the classical Kolmogorov one, we 
argue in Sect. 5 that it is not unrealistic on the parsec scale of 
the NGC 2264-C clump. 

Comparison between the observed dust continuum maps of 
NGC 2264-C and the synthetic column density maps of the 
"best-fit" simulation (see Fig. [6]) shows that the number of frag- 
ments and their alignment are well reproduced. By analogy 
with the observations, we have labelled the three main frag- 
ments of the synthetic 30m column density map SIM2, SIM3 
and SIM4. The corresponding synthetic PdBI map shows a 
strong central source, SIM3, surrounded by weaker sources, as 
observed. Moreover an additional component, labelled SIM13, 
becomes visible next to SIM3 in the simulations when "ob- 
served" at the PdBI angular resolution, which is strongly rem- 
iniscent of the (C-MM3, C-MM13) system in the real interfer- 
ometric map. 

As described in § 3.2, the collapse of the filament proceeds 
in two main phases: first, a global contraction velocity field is 
established along the long axis; second, a strong shock is gen- 
erated at the center by the two interacting sides of the model 
clump. At the same time, the clump fragments to form proto- 
stars. Thus, there are at least two relevant dynamical timescales 
in the problem: a global dynamical timescale corresponding to 
the global evolution of the elongated clump, and a local dy- 
namical timescale corresponding to the dynamical evolution of 
individual fragments. In other words, there is competition be- 
tween local collapse (i.e. fragmentation) and global collapse. 
In our simulations, this is controlled by the ratio of thermal 



to gravitational energy, oq., and thus by n c and ro (whose val- 
ues are 2000 cirT 3 and 0.7 pc, respectively) since the kinetic 
temperature and the mass of the model clump are fixed (cf 
§ 3.1). If n c or ro are too small (i.e. the density structure ap- 
proaches rtH 2 K r~ 2 ; cf Eq.©), the individual fragments do not 
have enough time to collapse on their own before entering the 
central shock. Therefore only one, massive central core forms. 
Conversely, if n c and ro are too large, many protostars (and 
eventually stars) form before any significant large-scale veloc- 
ity field is established along the clump long axis. Furthermore, 
in the latter case, the number of fragments produced in the 
simulations becomes larger than the observed number of frag- 
ments. 

Note that the collapse simulations shown in this paper were 
performed with model clumps of total mass M to , = 1000M o , 
while the total mass of NGC 2264-C is estimated to be some- 
what larger, ~ 1600 M Q (cf. PAB06). With more massive model 
clumps, we did not manage to reproduce the overall morphol- 
ogy of NGC 2264-C, in the sense that too much fragmentation 
occurred. This suggests that our models lack some source of 
support against gravity compared to the actual NGC 2264-C 
clump. This will be discuss further in § 5. 

4.2. The sharp central velocity discontinuity 

One of the most interesting features of NGC 2264-C is the cen- 
tral velocity discontinuity observed by PAB06 in optically thin 
tracers toward C-MM3 (see Fig.[7]l. This velocity discontinuity 
is believed to trace the axial (i.e., ID) collapse of NGC 2264-C 
along its long axis, as well as a possible dynamical interac- 
tion between protostellar sources at the center of the clump. 
Our new PdBI observations, which confirm the presence of a 
strong velocity gradient along the long axis of the clump (i.e. 
~ 8.4 km.s^'.pc -1 - see Fig. [H, set additional constraints on 
the velocity field in the central part of NGC 2264-C. 

When the initial level of turbulence was lower than 20%, 
our SPH simulations convolved to the 30m resolution led to 
a central discontinuity resembling that observed. Furthermore, 
the shape of the PV diagram observed at the PdBI resolution 
along the long axis of the clump appears to be a key tracer of 
the time evolution, as can be seen in Fig. [7] In addition to ob- 
served PV diagrams (first row), Fig. [7] shows the synthetic PV 
diagrams of our "best-fit" simulation (a!' , = 5%) convolved 
to the 30m resolution (first column) and to the PdBI resolution 
(second column) at the same four time steps as in Fig.[5](rows 
two to five). Note that simulations adopting a Kolmogorov-like 
turbulent energy spectrum and low values of or?. (< 20%) led 
to velocity discontinuities which are similar to the discontinu- 
ity shown here for the "best-fit" simulation. 

The reference time in Fig.|7]is the same as in Fig.|3J namely 
one global free-fall time (tff = 9.5 x 10 5 yr) after the start of 
simulation. At the first time step shown, tyy - 4 x 1 5 yr (second 
row), no clear kinematical signature is apparent in the synthetic 
PV diagrams, either at the 30m or at the PdBI resolution. At 
tff - 2 x 10 5 yr (third row), the synthetic PV diagrams start to 
exhibit a velocity structure reminiscent of the observed velocity 
gradient and central discontinuity, but the amplitude of the ve- 
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locity structure is not large enough to match the observations. 
At the best-fit time step, i.e., trt - 1 x 10 5 yr, the agreement 
between the simulated PV diagrams (fourth row) and the ob- 
served PV diagrams (first row) is quite remarkable. The central 
amplitude (i.e. ~ 2 km.s~'), shape, and position of the veloc- 
ity discontinuity are well reproduced. Moreover, the synthetic 
PdBI PV diagram shows a ~ 1 km.s~' velocity gap between the 
two central fragments, SIM3 and SIM13, as observed between 
C-MM3and C-MM13. 

The fifth row of Fig.|7]shows a later time, i.e., t// - 0.4 x 
10 5 yr, when the central shock is well developed. While the 
synthetic 30m PV diagram remains satisfactory, the synthetic 
PdBI PV diagram differs markedly from the observations. 
Note, however, that this late phase of evolution may not be cor- 
rectly described by our numerical model. Indeed, our simula- 
tions do not include feedback from protostars which clearly in- 
fluences the late dynamical evolution of cluster-forming clouds 
(cf. Li & Nakamura 2006). Thus, it is not clear if, in reality, 
the central shock would have time to develop as much as in the 
simulations to produce a PV diagram such as the one shown at 
the bottom right of Fig. [7] 

The time evolution of the a tur b ratio in the "best-fit" simu- 
lation is also given on the right-hand side of Fig. [7] While the 
initial level of turbulence was only or? , = 5% in this simu- 
lation, it can be seen the ratio of nonthermal kinetic energy to 
gravitational energy quickly increases up to a tur b = 27% at the 
"best-fit" time step and a, ur b = 33% at the last time step shown. 
This demonstrates that the bulk of the "turbulent" energy in our 
simulations does not come from the large scale turbulent veloc- 
ity field, but rather from the conversion of gravitational energy 
into kinetic energy through the global collapse of the clump. 
The increase of a tur b with time contributes to produce synthetic 
linewidths in reasonable agreement with observed linewidths 
(see § 4.3 and Fig. [8]below), despite the low level of kinetic 
energy at the beginning of the simulation. We also note that 
the value of a tur b achieved at the "best-fit" time step (27%) is 
within a factor of 2 of the kinetic to gravitational energy ra- 
tio expected in virial equilibrium (50%), despite the fact that 
the model clump is globally collapsing and far from equilib- 
rium at this stage. Clearly, the broad linewidths observed in 
NGC 2264-C result at least partly from systematic inward mo- 
tions as opposed to random turbulence. In our "best-fit" model, 
most of the motions are gravitationally focussed and do not ex- 
ert any support against gravitational collapse. 

4.3. Large-scale kinematical pattern in optically-thin 
line tracers 

The H 13 CO + (1-0) spectra observed toward NGC 2264-C with 
the 30m telescope show a remarkable East- West axial symme- 
try over the whole clump on either side of the central source 
C-MM3 (see Fig.©. The low optical depth of the H I3 CO + (l- 
0) line (r ~ 0.3 for the peak velocity channel of the central 
spectrum) inferred from the Monte-Carlo radiative transfer cal- 
culations of PAB06 implies that the observed double-peaked 
line profiles (cf. Fig. [8} result from the presence of two veloc- 
ity components along the line of sight rather than from self- 



absorption. Although our PdBI observations show that the cen- 
tral velocity discontinuity seen with the 30m telescope origi- 
nates from two unresolved protostellar sources, the large spatial 
extent of the region over which double-peaked H I3 CO + (1-0) 
spectra are observed in the 30m map suggests a global kine- 
matical origin for the double-peak profiles. 

While the N 2 H + (101-012) PV diagram observed with PdBI 
(Fig.|4]i may be suggestive of rotation about an axis perpendic- 
ular to the long axis of the clump, PAB06 showed that rotation 
could not account for the shape of the observed 30m PV di- 
agram based on a detailed comparison with radiative transfer 
models (see Fig. 12 of PAB06). By contrast, we now proceed to 
show that our scenario of large-scale, axial collapse does pro- 
vides a good match to the symmetric pattern of double-peaked 
line profiles observed in low optical depth tracers. 

To this aim, synthetic spectra were constructed from our 
SPH simulations assuming strictly optically thin line tracers: 
each SPH particle was given the same weight and the contri- 
butions of all particles falling within a given velocity channel 
were integrated. The synthetic data cube was then convolved to 
the 30m angular resolution and normalized in such a way that 
the peak intensity of the synthetic central spectrum matched 
the peak intensity of the observed central H I3 CO + (1-0) spec- 
trum. Figure |8]compares the H I3 CO + (1-0) spectra observed in 
the central part of NGC 2264-C (top) with the resulting syn- 
thetic spectra for the best-fit simulation at the best-fit time step 
(bottom). (In this comparison, we use the observed H 13 CO + (l- 
0) spectra rather than the observed N2H + (l-0) spectra since the 
former have a better signal-to-noise ratio.) It can be seen that 
the overall agreement is very good. Since the synthetic line 
emission is optically thin, the double-peaked spectra exhibited 
by the model are clearly not due to radiative transfer effects but 
result from the presence of two velocity components along the 
line of sight, corresponding to the two ends of the elongated 
clump moving toward each other. Focusing on the central row 
of spectra (marked in boldface), it can be seen that the blue- 
shifted component of the double-peaked spectra dominates on 
the eastern side. Moving west, the red-shifted component be- 
comes progressively stronger. It is nearly as intense as the blue- 
shifted component at the central position and eventually dom- 
inates on the western side of the filament. This remarkable re- 
versal of blue/red spectral asymmetry as one moves from the 
eastern to the western side of the central C-MM3 position can 
be seen in both the observations and simulations. 

The synthetic spectra obtained from our SPH simulations 
are mass weighted and are thus more representative of the 
global kinematics of the clump than of the kinematics of com- 
pact individual fragments. We conclude that the remarkable 
pattern seen in the central row of spectra in Fig. [8] character- 
izes the collapse of the elongated clump along its long axis. 
We note, however, that the synthetic spectra are somewhat nar- 
rower than are the observed H 13 CO + (1-0) spectra. Part of the 
observed linewidths may result from outflowing gas generated 
by the protostars, an effect which we did not treat in our simu- 
lations. It may also be partly due to another source of support 
against gravity, not included in our simulations (see § 5). 

In the context of our interpretation of the double-peaked 
spectral pattern observed in H 13 CO + (1-0), the extent of the re- 
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gion over which the red-shifted peak is observed on the eastern 
side (and the blue-shifted peak observed on the western side) 
sets constraints on the diameter of the NGC 2264-C cylinder- 
like clump (see Fig. 1 1 of PAB06). Double-peaked H 13 CO + (l- 
0) spectra are observed up to 30" on either side of the cen- 
tral object C-MM3. Given the distance of 800 pc and assuming 
a viewing angle of 45 degrees between the line of sight and 
the long axis of the clump (as adopted in the radiative transfer 
model presented by PAB06), we estimate the diameter of the 
cylinder to be ~ 0.65 pc. This is in good agreement with the 
apparent width of the NGC 2264-C clump as measured in the 
plane of sky on our dust continuum and molecular line maps. 

5. Concluding remarks 

The good quantitative agreement obtained between our "best- 
fit" SPH simulations and our (30m and PdBI) millimeter ob- 
servations confirms the physical plausibility of the scenario 
of large-scale axial collapse and fragmentation proposed by 
PAB06 for the NGC 2264-C clump. Observationally, such an 
axial collapse is traced by a central velocity discontinuity as- 
sociated with double-peaked profiles in optically thin line trac- 
ers. The present study supports our earlier suggestion that an 
ultra-dense protostellar core of mass up to ~ 90 M Q is in the 
process of forming at the center of NGC 2264-C through the 
dynamical merging of lower-mass Class cores (cf. PAB06). 
Our interferometric PdB detection of a new object, C-MM13, 
located only ~ 10000 AU away (in projection) from the cen- 
tral source, C-MM3, but with a line-of-sight velocity differing 
by ~ 1.1 km.s~' from that of C-MM3, provides an additional 
observational manifestation of the merging process. Given the 
relatively large mass of C-MM13 (~ 8M ), such a large ve- 
locity difference would be difficult to explain by dynamical 
fragmentation during the collapse of an individual protostellar 
core, even if low-mass objects can easily be ejected from dy- 
namically unstable protostellar systems (e.g. Bate et al. 2003, 
Goodwin et al. 2004). In our proposed scenario for NGC 2264- 
C, the local collapse of individual protostellar cores is strongly 
influenced by the high dynamical pressure resulting from the 
global collapse of the clump, and proceeds in a manner that is 
qualitatively similar to the triggered protostellar collapse mod- 
els discussed by Hennebelle et al. (2003, 2004). 

Our detailed comparison between observations and simula- 
tions has also allowed us to set constraints on the evolutionary 
state of the NGC 2264-C protocluster. It seems that the charac- 
teristic shape of the observed position-velocity diagrams sur- 
vives for a relatively short period of time, i.e. < 1 x 10 5 yr, and 
occurs only very soon after the formation of the protocluser 
while less than ~ 1% of the gas has been accreted onto sink 
particles. 

The low level of initial turbulent energy required to match 
the observations implies that NGC 2264-C is structured more 
by self-gravity than by turbulence. The main effect of turbu- 
lence is to create seeds for further gravitational fragmentation. 
Turbulent fragmentation does not appear to play a significant 
role in this clump. In our "best-fit" simulation, the initial tur- 
bulent to gravitational energy ratio is a Q turb = 5%, comparable 
to the ratio of thermal to gravitational energy aP,. The level of 



"turbulence" increases as the simulation proceeds and gravita- 
tional energy is converted into kinetic energy. At the "best-fit" 
time step, the ratio of nonthermal kinetic energy to gravita- 
tional energy approaches ~ 30%. Most of the corresponding 
"turbulence" is gravitationally generated as in the recent cloud 
collapse simulations of Burkert & Hartmann (2006). In other 
words, the cloud motions in our best-fit model are primarily due 
to collapse and gravitationally organized motions as opposed to 
purely random turbulence. Although we have not identified a 
specific trigger, we believe that the "cold" or "subvirial" initial 
conditions (cf. Adams et al. 2006) required by our model reflect 
the fact that the NGC 2264-C clump was suddenly compressed 
and/or assembled as a result of a strong external perturbation. 

The fact that simulations starting from initial turbulent ve- 
locity fields with a Kolmogorov-like energy spectrum lead to 
model clumps that are much less organized than the observa- 
tions should not be overinterpreted. Indeed, since the phases 
are chosen randomly and since in Kolmogorov-like turbulence 
most of the energy is on large scales, it is not surprising that the 
shape of the filament is strongly distorted in this case. In a real 
situation, the large-scale turbulent fluctuations should be much 
more coherent since they may be responsible, at least in part, 
for the formation of the filament in the first place (cf. Hartmann 
et al. 2001). 

Another point worth noting is that the total mass of gas with 
density above 10 4 cirT 3 is ~ 10 times lower in our best-fit sim- 
ulation than in the actual NGC 2264-C clump. Using higher 
densities by a factor of 10 in the numerical simulations would 
inevitably lead to fragmentation into a larger number of cores 
since the corresponding Jeans mass would be smaller by a fac- 
tor ~ 3 compared to the Jeans mass in the present simulations. 
It seems therefore that some additional support against gravity, 
not included in the simulations presented here, plays a role in 
NGC 2264-C. This extra support could arise from protostellar 
feedback and/or magnetic fields. 

Finally, we speculate that the evolution inferred and simu- 
lated here for NGC 2264-C is not exceptional but representative 
of many massive cluster-forming clumps in the Galaxy. In par- 
ticular, we note that evidence of large-scale, supersonic inward 
motions has been recently found in several deeply embedded 
regions of high-mass star formation (Motte et al. 2005 - see 
also Wu & Evans 2003 and Fuller et al. 2005). NGC 2264-C 
may just be caught at a particularly early stage of protocluster 
evolution and observed in a favorable configuration, leading to 
a remarkably simple kinematical pattern. Similar detailed mod- 
elling studies of other cluster-forming clumps will be needed to 
confirm this hypothesis. 
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Fig. 6. Observed column density distribution (first row) compared to synthetic column density maps (second to fifth row) con- 
volved to the 30m angular resolution (first column) and to the PdBI angular resolution (second column) for four different initial 
levels of turbulence: or b — 1% (second row), a° turb = 5% (third row), a® b = 20% (fourth row), aft , = 50% (fifth row). The 
best-fit simulation corresponds to the third row, i.e., aP b = 5% (for which the displayed time step is the "best-fit" time step). 
Note that the synthetic PdBI maps include the effect of interferometric filtering (see text). In each map, the contour levels go 
from 10 to 90% by step of 10% of the peak emission. The observed 30m column density distribution corresponds to the 1.2mm 
dust continuum map of PAB06. 
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Fig. 7. Observed position-velocity diagrams (first row) compared to synthetic position-velocity diagrams (second to fifth row) 
convolved to the 30m angular resolution (first column) and to the PdBI angular resolution (second column) at four different time 
steps of our best-fit SPH simulation (a Q turb = 5%): t// - 400000 yr (second row); t// - 200000 yr (third row); t// - 100000 yr 
(fourth row); t// - 40000 yr (fifth row) (with t// = 9.5 x 10 5 yr). The best-fit time step corresponds to fourth row. The value of 
the turbulent to gravitational energy ratio a tur b at each time step is given on the right hand side. Note that a, m -b increases as time 
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NGC 2264-C: 30m observations 
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Fig. 8. Comparison between the H 13 CO + (1-0) spectra observed at the 30m telescope (upper panel) and the synthetic optically 
thin spectra obtained in our best-fit simulation (lower panel). In the upper panel, the (0,0) corresponds to the position of C-MM3, 
while in the lower panel the (0,-10) position corresponds to the position of SIM3. Overlaid in grey scale are the 1.2 mm dust 
continuum image of PAB06 (top) and the synthetic column density map of the best-fit simulation (bottom). The row of spectra 
observed (top) and simulated (bottom) along the main axis of the clump are marked in boldface. 



